

** make bootstrap cis
/*
1. append together bootstrap sample
2. in each sample-week, calcualte bounds based on non-icli and clear cause hospitalizations
3. find 2.5th and 97.5th percentile of bounds (by week) 
4. find average lb and ub, by week, and get weekly bias in bound 
*/

** 1. append together
local B 500
clear
forvalues bs = 0/`B' {
	append using bootstrap_ests/bs_`bs'
}
save bootstrap_ests/bs_all, replace

** 2. calculate intersection bounds , ub_non_m only
use bootstrap_ests/bs_all, clear
gen ub_non_m = min(upper_bound1, upper_bound2)
gen ub_clear_m = min(upper_bound1, upper_bound3)
gen ub_icli_m = min(upper_bound1, upper_bound4)

sum week upper_bound2 if bs == 0

local n 0
local weeks `=mdy(4,24,2020)' `=mdy(6,19,2020)' `=mdy(3,27,2020)'
foreach week of local weeks {
	local ++n
	sum ub_non_m if week==`week' & bs == 0
	local estimate = string(r(mean), "%3.2f") 
	
	if `week' == `=mdy(4,24,2020)'  local name April 24
	if `week' == `=mdy(6,19,2020)'  local name June 19
	if `week' == `=mdy(3,27,2020)'  local name March 27
	
	kdensity ub_non_m if week==`week' & bs ~= 0, nograph gen(x d) 
	sum x
	local xp = r(min)
	sum d
	local yp = r(max)
	# delimit ;
	twoway
		(line d x)
		, 
		xtitle("Estimate") ytitle("") 
		graphregion(color(white))
		text(`yp' `xp' "Week of `name'" "Estimate: `estimate'", 
			place(e) justification(left) )
		name(g`n', replace)
	;
	# delimit cr
	drop d x
	
	
}

graph combine g1 g2 g3, cols(1) graphregion(color(white)) xsize(4) ysize(9)
graph export figures/density_bs.pdf, replace 
